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We study discrete solvent effects on the interaction of two parallel charged surfaces in ionic 
aqueous solution. These effects are taken into account by adding a bilinear non-local term to the 
free energy of Poisson-Boltzmann theory. We study numerically the density profile of ions between 
the two plates, and the resulting inter-plate pressure. At large plate separations the two plates are 
decoupled and the ion distribution can be characterized by an effective Poisson-Boltzmann charge 
that is smaller than the nominal charge. The pressure is thus reduced relative to Poisson-Boltzmann 
predictions. At plate separations below ~ 2nm the pressure is modified considerably, due to the 
solvent mediated short-range attraction between ions in the the system. For high surface charges 
this contribution can overcome the mean-field repulsion giving rise to a net attraction between the 
plates. 



I. INTRODUCTION 

Aqueous ionic solutions are abundant in biological and 
chemical systems. Often they play a prominent role in de- 
termining the properties of charged macromolecules that 
are immersed in thenJj. The mean field theory of elec- 
trolytes, known as Poisson-Boltzmann (PBl. theory and 
its linearized version, Debye-Hiickel theorytl~l3, are known 
for many decades and have proved to be useful and im- 
portant tools. PB i^keory was applied in tbj3 study of 
colloidal dispersionaD'El, biologiC|al.membranesEl, synthetic 
and biological polyelectrolytesElJ, and complex systems 
such as DNA-lipid complexesEJ. Nevertheless, PB the- 
ory is known to have important limitations. Being a 
mean field theory, ion-ion correlations are ignored. In 
addition, the finite size of ions is neglected. These effecis. 
have been studied e»teiiS|ively using various apprpachestj 
such as liquid stotalp~ll3 and densiiy.. functional theo- 
ries, simulationaiaEj, fieid tiieoryEilH and other modifi- 
cations to the PB theoryl23~E3. 

Most of the studies of corrections to PB have con- 
centrated on the so-called primitive model, where ions 
are assumed to interact with each other through the 
electrostatic interaction and a hard core steric repul- 
sion. Although this model can describe many effects 
that are neglected in PB theory, it still neglects some 
physical features that are present in real systems. Most 
notably, the aqueous solvent is treated as a continuous 
medium, whereas in reality ions interact with discrete 
solvent molecules. 

Solvent effects are strong especially in water, because 
the polar water molecules interact very strongly with 



ions. The most significant result is that the electrostatic 
ion-ion interaction is reduced by a factor e ~ 78 at room 
temperature, due to screening by the dielectric environ- 
ment. However the discreteness of the solvent results in 
a more complicated picture. When ions approach each 
other at separations of a few water molecular diameters, 
the effective interaction between them is modified consid- 
erably. Fig. 1 shows the correction to the 1/er potential 
between two Na+ ions in water. This effective potential 
was calculated, using a simulation schemecZI, for a bulk 
NaCl solution of concentration 0.55 M, at room temper- 
ature. Note that the short-range potential, remaining 
after the subtraction of the Coulomb interaction, is os- 
cillatory and predominantly attractive. 

The possibility to calculate the effective potential be- 
tween ions in water leads naturally to the model depicted 
schematically in Fig. 2. The water is treated as a contin- 
uous medium, with a dielectric constant e. In addition 
to the electrostatic interaction a short-range interaction 
is included between ion pairs. The short-range potential, 
denoted (r), is taken as an input to the model (from 
simulation), and can in general depend on the ion species 
i and j. For example, the potential shown in Fig. 1 is used 
between Na'^'-Na"'" pairs. The effective potential is cal- 
culated in a bulk solution and thus depends only on the 
ion-ion separation. However, systems containing charged 
surfaces can lead to inhomogeneity or anisotropy in the 
ion distribution. 

The model described above was suggested in 
Ref. |2^, and was studied in planar geometry us- 
ing the Anisotropic Hyper-Netted Chain (AHNC) 
approximatiorO in Refs. In Ref. p3 we presented 
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a simplified approach to the same model In this latter 
approach, a term accounting for the short-range solvent- 
mediated ion-ion interaction is added to the PB free en- 
ergy. The formalism obtained in this way is simple al- 
though less accurate than the AHNC approximation, and 
in particular neglects ion-ion correlations. On the other 
hand numerical calculations can be done fairly easily, and 
are feasible in non-planar geometries. In addition, vari- 
ous analytical results can be obtained, and the discrete 
solvent effects can be readily understood in terms of ba- 
sic physical principles. In the present paper, which can 
be regarded as a follow-up of Ref. |3^, we use the same 
formalism to study discrete solvent effects on interacting 
charged and planar plates. 

The outline of the paper is as follows. Section || 
reviews the model and discusses its application to two 
charged and planar plates. In Sec. Ill we discuss the cor- 
rections to the PB density profile. In section IV we obtain 
expressions for the inter-plate pressure and derive a gen- 
eralized contact theorem. The resulting pressure curves 
are studied numerically and analytically in section 
Finally, section VI offers some concluding remarks. The 
technical details in the derivation of the pressure are pre- 
sented in the Appendix. 



The first three terms in Eq. form the usual PB 
expression for the free energy. The first term is the elec- 
trostatic free energy and the second term is the entropy 
of the ions. The electrostatic potential 5" is a functional 
of the ion densities Ci, and is determined by the Poisson 
equation and the boundary conditions imposed by the 
surface charges. Instead of writing this dependence ex- 
plicitly in the free energy, it is convenient to add a third 
term to 57, containing a Lagrange multiplier A(r). 

The fourth term in Eq. (|l|) accounts for the hydration 
interaction, and is quadratic in the ion densities. The 
weighted potential Uij is defined as: 



U,, = 1 



-/3«,,(|r-r'|) 



(2) 



where Uij is the nominal short-range hydration interac- 
tion between ions of species i and j. To obtain Eq. (|l|) 
we first treat the short-range interaction Uij using a virial 
expansion of the grand canonical potential, keeping terms 
up to the quadratic order. The electrostatic interaction 
is then treated exactly as in PB theory, using a mean 
field approximation for the electrostatic potential As 
an alternative approach Eq. (0) can be obtained frora-a 
field theory expansion of the grand partition functions. 



II. THE MODEL 

A. Free energy 

The free energy of the system can be written as a func- 
tional of the local ion densities, consisting of the usual 
PB term and a hydration correction term. Assuming 
that the boundary conditions are of fixed charges, the 
following-approximated form for the free energy can be 
obtained^: 

n=^l {V^)M'r + kBT jY. (\n ^ - d^r 

+ j k{v) (^H + ^Y.^,e}j d^r 

+ / c.(r)c,(r')f/.,(r-r')d3rd3r' (1) 

where ^ is the electrostatic potential, Ci are the ion den- 
sities, Ci are their respective charges, e is the dielectric 
constant, ksT is the thermal energy and Uij is defined 
below. The bulk ion densities Cb^i are determined by the 
fugacities C,i — exp(/3//i)/Ai^, where jii are the chemical 
potentials, is the de Broglie thermal wave length and 
(3 = l/fcsT. The simple PB relation Cb_j = Qi is altered 
with the inclusion of hydration interactions, as will be 
explained below. A detailed discussion of the various ap- 
proximations involved in Eq. is given in Ref. Here 
we shall briefiy discuss each of the terms, and outline the 
way in which Eq. (|^) is obtained. 



B. Density equations 

The density profiles are obtained by minimizing the 
free energy 57 with respect to the ion densities c^. The 
third term in Eq. (^, containing the Lagrange multiplier 
A(r) allows us to regard the densities Ci(r) and the elec- 
trostatic potential ^'(r) as independent fields, and require 
that 57 has an extremum with respect to the three fields 
Ci, ^E" and A. Requiring that 57 has an extremum with 
respect to ^' gives: 

A^i* (3) 

and the extremum condition with respect to Ci then gives: 
In + ^ /■ c,(r')C/., (r - r') d^r' + /3e,vl/(r) = 

(,^ ^ J 

(4) 

where relation (^) has been substituted to express A in 
terms of 5*. This equation is supplemented by the Pois- 
son equation: 

VH = -^Y.^,c. (5) 

i 

Since Eq. (^) is an integral equation, the Ci cannot be 
written as a simple function of Therefore a single 
equation for 5*, analogous to the PB equation, cannot be 
obtained, and we are left with the two coupled equations 
(Q) and (H). These equations should be solved together 
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to obtain the electrostatic potential and density profiles. 
For U ^ 0, Eq. (0) reduces to the Boltzmann equation 
Ci — Ci exp(— /3ei^. In the bulk = 0, leading to the 
relation Ch^i — Combining these relations with Eq. (j^) 
reproduces the PB equation: 



£ ^ 



(PB) 



(6) 



Equations (||) and (||) were solved for a single charged 
and planar plate in Ref. The treatment of two parallel 
plates is very similar, and is outlined below for complete- 
ness. The system is shown schematically in Fig. 2b. The 
plate positions are designated by 2; = and z ~ d, using 
the convention that these are the coordinates of closest 
approach of the ions to the plates (while the potentials 
Uij(r) are measured from the centers of the ions). The 
two plates are negatively charged, each one with a uni- 
form surface charge a. No discreteness of surface charge 
is taken into account in the present work. We assume an 
electrolyte of valency i.e., a solution of positive 

and negative ions of charges e± — ±z±e, where e is the 
electron charge. 

In order to simplify the equations further, the inter- 
actions between the different pairs of ion species can be 
taken to be equal, i.e., Uij{r) = U++{r) = U{r) where 
is the weighted potential between the (positive) 
counterions. The exact choice of and U^^ is ex- 
pected to be of only minor significance, as the co-ions are 
repelled from the surface neighborhood and only the pos- 
itive counterions reach high densities there. From charge 
neutrality we have Cf, = Cf,_+ = (z_/z+)cf,^_ and similarly 

= = (z_/z_|_)C_, where the relation between Cb and 
C will be determined later. 

Due to the one-dimensional symmetry imposed by the 
charged and planar planes, the integration in Eq. (|j) can 
be performed over the x — y plane to obtain: 



c±{z) = C±e 



exp 



c{z')B(z — z') dz' 



(7) 



where c = c+ -|-c_ is the total ion density and B(z) is the 
effective interaction between two layers of ions, expressed 
as an integral of U{r) in the plane of constant z. Using 
cylindrical coordinates: 



/•OC 

B{z) = 2tt pdpUi^/z^ + p^) 
Jo 



(8) 



The Poisson equation reads: 
d^\P 47re 



dz2 



— (e 



X exp 



c{z')B{z - z')dz' 



(9) 



Equations (^ and (|^) are supplemented by the following 
boundary conditions: 



d* 
d7 



An d^P 
= cr ; — 

£ dz 



= 



(10) 



z=d/2 



since the problem with two plates of equal charge at z = 
and 2 = d of equal charge is symmetric about the mid- 
plane z — d/2. 

Finally, the relation between ^ and the bulk density 
Cb can be obtained from Eq. (^. We imagine that the 
two plates are immersed in a bath of electrolyte. In the 
region outside the plates an equation similar to Eq. (|^) 
holds, where the integration inside the exponent is per- 
formed in the external region. Far away from the plates, 
as 5* becomes zero, c+ and c_ assume their asymptotic 
constant, bulk values. The integrand inside the exponen- 
tial can be replaced by —(1 -I- z+/z-)cbB{z — z') leading 
to the result: 



Cb = C exp 



where: 



Bt= f dzB{z) = I d^YU{v) 

J —00 J 



(11) 



(12) 



is also equal to 2i?2, the second virial coefficient. The 
limit BtCb ^ is the limit in which the short-range in- 
teraction becomes negligible in the bulk. In this limit the 
relation between the bulk density and fugacity of Eq. ( pT| ) 
tends to the ideal gas relation Cb = C = exp(/3^)/Aj.. 

In the next section we will concentrate on a symmetric 
1 : 1 electrolyte, where the equations (|^) and @ take the 
form: 



c±(^) = Ce^'^"*exp 
d^^f 8ne 



c{z')B{z - z')dz' 



dz2 



Csinh (/3e^') 

£ 



X exp 



c(z')B{z — z') dz' 



and: 



Cb = Cexp {-2BtCb) 



(13) 



(14) 



C. Definitions and parameters 

For the short-range ion-ion potential u{y — r') we 
use the effective potential between Na^ - Na"^ ion pairs, 
shown in Fig. 1. For ion-ion separations below 2.9A a 
hard core interaction is assumed. Fig. 3 shows the ef- 
fective layer-layer interaction B{z), as was derived from 
this potential using Eq. (||). This effective interaction 
is mostly attractive, as B{z) is negative on most of its 
range, and has a characteristic range of approximately 



3 



7 A. The structure of B{z) reflects the oscillatory behav- 
ior of u{r). 

It is useful to introduceyJ^he length scales characteriz- 
ing the PB density profileaj. The Gouy- Chapman length, 
defined as 6 = eksT / {2TTe\a\), characterizes the width 
of the diffusive counterion layer close to a single charged 
plate with a surface charge density cr, in the absence of 
added salt. The Debye-Hiickel screening length, Ad ~ 
{SiTCbe^/ekBT)-^/^, equal to 19.6 A for Cb = 0.025 M at 
room temperature characterizes the decay of the screened 
electrostatic interaction in a solution with added salt. 
The strength of the electrostatic interaction can also be 
expressed using the Bjerrum length, Zb = /{eksT). 
This is the distance at which the electrostatic interaction 
between two unit charges in a dielectric medium becomes 
equal to the thermal energy. It is equal to about 7A in 
water at room temperature. In terms of the Bjerrum 
length h = e/27rZB|(T| and \d — (SncblB)'^^^ ■ 

The inclusion of the hydration interaction introduces 
additional length scales in the system. For the interac- 
tion of Figs. 2 and 3, the range of the interaction dhyd 
is approximately 7 A, over twice the hard core diameter 
dhc = 2.9A. The strength of the hydration interaction is 
characterized by the second virial coefhcient Bt/2, with 
Bt ~ -(7.9A)3 as is calculated from Eq. (H). 



III. DENSITY PROFILES 

Equations (0) and (||) are a set of three nonlinear in- 
tegrodifferential equations. We treat them numerically 
using an iterative scheme, based on the assumption that 
the positive ion density profile is dominated by the elec- 
trostatic interaction. We start with the PB profile and 
calculate iteratively corrections to this profile, as result 
from Eqs. (0) and (|). For a 1:1 electrolyte we iteratively 
solve the equation: 



d2*(") 8 



dz2 



^Csinh(/3e*(")) 



X exp 



=("-i)(z')S(2-z')d2' 



(15) 



where the superscript n stands for the nth iteration, 

4")(.)^Ce^ 



X exp 



"J''-^\z')B{z - z')dz' 



(16) 



and the zeroth order densities Cj? ^ are taken as the density 
profiles generated by the PB equation (^) . The boundary 
conditions ( |To| ) are satisfied by the electrostatic potential 
in all iterations. The solution converges after several it- 
erations. It is interesting to note that the first iteration 
captures most of the effect. This observation can lead to 



various analytical results, as shown in Ref. 33 for a single 
plate. 

In the following sections we will concentrate on the 
pressure between the plates. First we discuss briefly the 
modification to the PB density profiles. Let us begin by 
considering a large plate separation d. In this case the re- 
sults are similar to the single-plate case, since d is larger 
than all other length scales in the system, and we present 
them for completeness. 

Fig. 4 shows the density profile of the positively 
charged counterions (solid line) between two charged 
plates, with d = 50A. Only one half of the system is 
shown, since the profile is symmetric around the mid- 
plane. The surface charge, \a\ = 0.333 C/m^ corresponds 

to approximately 48 A per unit charge. This is a typical 
high surface charge obtained for mica plates. It corre- 
sponds to a Gouy-Chapman length b = 1.06 A, at a tem- 
perature of 298 K, with e = 78. The electrolyte bulk con- 
centration is 0.025 M, corresponding to a Debye-Hiickel 
screening length Ad = 19.58 A . The density profile is 
compared to the result of PB theory (dotted line). 

The main effect is that the short-range attraction 
draws additional counterions to the vicinity of the 
charged plate. Note, however, that the contact density 
remains very close to the PB density, as will be explained 
later. The increase of the counterion density near the 
plate is followed by a depletion further away. This can 
be understood in the no-salt case since the total number 
of counterions is fixed. In our case the salt concentra- 
tion is low. The Debye-Hiickel screening length is large 
compared to the Gouy-Chapman length and compared 
to the range of the short-range interaction, so the salt 
has a minor effect. 

The counterion density profile ia_also compared with 
results of the AHNC approximationEj that were obtained 
using the same short-range hydration potential ('x' sym- 
bols). The qualitative effect is similar in our model and 
in the AHNC. Specifically, both density profiles follow 
the PB density curve for the first few Angstroms from 
the plate and show a considerable decrease in the posi- 
tive ion density, relative to PB, starting at a distance of 
about 5 A from the plate. The maximal decrease in the 
density is approximately 30% in our model and almost 
50% in the AHNC profile, both relative to the PB profile. 

The effect of the short-range ion-ion interaction 
strongly depends on the surface charge. This is demon- 
strated in Fig. 5. The ratio of the counterion density and 
its PB value, c/c^^, is shown for three values of a. The 
effect of the hydration potential is very minor for small 
surface charge (|cr| = 0.0333 C/m^ ~ le/480A^), where 
the ratio c/c^^ is approximately 2% at its maximum, 
and considerable for a surface charge of 0.333 C/m-^ = 

1 e/48 A , where it reaches approximately 40%. 

As the plate separation decreases, the modification to 
is expected to remain similar to the single plate case 
as long as d/2 is large compared to b and to c?hyd- This 
can indeed be seen in Fig. 6, where a high surface charge. 



„PB 
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as in Fig. 4, is considered. In this case b ^ dhyd, so a 
deviation from the single plate curve is expected when 
d/2 ^ dhyd — 7 A. The ratio c/cpB is shown for several 
plate separations between 5 and 50 A. The results are 
very similar for d = 50, 35 and 20 A (Fig. 6a). In partic- 
ular, note that the contact density remains very close to 
the PB value in all three separations. This is a result of 
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the generalized contact theorem, derived in section IV. 
For smaller d, 5 and 10 A (Figures 6b and 6c, respec- 
tively) the behavior is different, and in particular the 
contact density deviates from the PB value. The effect 
of decreasing d was found to be similar for smaller surface 
charge {e.g., 0.1 C/m^,) and for salt concentration up to 
0.1 M. 

The most important effect on the density profile is that 
the ion density is depleted far away from the charged 
plates. When the two plates are highly separated from 
each other the ion density can be described, far away from 
the plates, using an effective PB surface charge. This ef- 
fective charge was calculated in Ref. 3^, and is smaller 
than the nominal charge (for example, for the surface 
charge used in Fig. 4 it is smaller by a factor of ~ 3.8). 
The reduced density leads to a reduced pressure, relative 
to PB, as will be explained in the following sections. 



IV. PRESSURE EQUATION AND CONTACT 
THEOREM 



The pressure Pi„ in the region between the two plates 
can be obtained by differentiating the free energy Q with 
respect to the plate separation d: 



Pr 



6n 

Id 



(17) 



To compute 6fl we can imagine that a 'slice' of width Sd 
is inserted at some position zq between the two plates. 
Adding up all the contributions to Sft, and using Eq. (|^) 
and the boundary conditions ( p^ we obtain: 



d 



kBTj^ dz dz'a,{z)cj{z') — {z'- 

Jo Jzo 



z) (18) 



This result is correct for any combination of ion species 
i, assuming the same short-range interaction Uij between 
different ion pairs. The full derivation is given in the 
Appendix. The pressure is equal throughout the plate 
spacing and, therefore, independent on the choice of zq. 

The net pressure P between the plates is the difference 
between the pressure inside and outside the plates. The 
latter is equal throughout the region outside the plates. 
In particular, it is equal to the bulk pressure Pbuik, so we 
have: 



bulk 



(19) 



To obtain Pbuik, we note that an equation similar to 
Eq. (|l^) holds in the bulk, with constant electrostatic 
potential and with Ci constant and being equal to the 
bulk densities. For the case of a 1:1 electrolyte, we find: 



Py 



bulk 



2kBTcb (1 + BtCb) 



(20) 



Since Bt is negative the bulk pressure is lower than its PB 
value. Note that in the case of no added salt Pbuik = 0. 



The expression (18) assumes a particularly simple form 
if we set zq to zero^namely on one of the plates. Then 
the third term in (|8|) vanishes and the second term is 
fixed by the boundary conditions, giving: 



p = fcsr^Q(o) 



27r 



Pi 



bulk 



(21) 



Alternatively, if we choose zq at the mid-plane, z = d/2, 
by symmetry the second term in ([l^) vanishes and the 
pressure is expressed as: 



P = kBTY,c^[d/2) 

fd/2 I'd 
Id/ 



/•d/2 pd 1 jo 

kBTj2j dz dz' c,{z)c,{z') — {z' - 



Phul] 



(22) 



The equality of these two expressions fop|j±iie pressure 
results in the generalized contact theoremlyE3: 



rd/2 pd 
I d/2 



Y^fJ^ dz j'^ dz' c, (z)c, [z')'^{z'-z) (23) 



The very small relative change of the contact density, 
compared to PB theory, at large plate separations can be 
understood from this result. We consider first the case of 
high surface charge, where the the Gouy-Chapman length 
is small compared to the Debye-Hiickel screening length, 
6 ^ Ad. In this case the second and third terms on the 



right hand side of Eq. (23) become negligible compared to 
the first term when d ^ b, dhyd, where dhyd is the range 
of the hydration interaction. The contact ion density is 
then dominated by the positive ion density, and is very 
close to the PB value. When there are only counterions 
in the solution and c? — s- oo (or equivalently, in the case 
of one isolated plate) , we have exactly, as in PB theory : 



c+(0) 



27r/3 2 
— -a^ 

e 



(one plate, no salt) 



(24) 



If b is not small compared to A^i, the correction to 
the contact density is still small for large enough plate 
separations, assuming that the hydration interaction is 
negligible in the bulk, i.e. —BtCh ~ — i?t/(87rZBAp) <C 1. 
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When d ^ X]j and d > dhyd, the couphng between the 
two plates is negligible and Eq. (123) becomes: 



Ci:(0) ^ —a 



bulk 



(25) 



The only difference in this expression relative to the PB 
contact density is the change in the bulk pressure. This 
change is negligible if the hydration interaction is negli- 
gible in the bulk. 

For smaller d, the integral in (^3|) can contribute to a 
significant change in the contact density relative to PB 
theory. This can be seen in Fig. 6, where b << Ad, at 
plate separations below 10 A. 



V. PRESSURE CURVES 



A. Pressure beyond Poisson-Boltzmann 



Using equations ( p2[ ) and (|2C|) the pressure can be writ- 
ten in our model as the sum of the following three terms: 

P = kBTJ2[c^{d/2)-Cb^,] 

i 

fd/2 rd 

iz 

Id/2 



* ' z — z) 



2kBTBtcl 



(26) 



A symmetric 1:1 electrolyte is assumed for simplicity 
throughout this section. We would like to compare this 
pressure with the PB pressure, which can be written as 
follows: 



PpB = fcsT^ [cpB,,;(d/2) - Cb,i] 



(27) 



where cpB,i(rf/2) is the PB density of the zth ion species 
at the mid-plane. The first term in Eq. (pF 



Pm = kBTj2h{d/2) - Cb,^] 



(28) 



is similar in form to the PB pressur e (27| ), but the mid- 
plane density in equations (|28| ) and (^) can be different. 
The second term in Eq. (Eq), which we denote as the 
hydration pressure: 



P 



hyd 



d/2 



dz' Ci{z)cj{z')— — (z' — z) 



dz 



(29) 



is the integrated short-range force acting between ion 
pairs in the two halves of the system. The third term 
is the change in the bulk pressure relative to PB theory, 
due to the inclusion of a 2nd virial coefficient in the bulk 
equation of state. 



Some simple observations can be made immediately 
from Eq. (^6|). These observations will be useful in the 
next subsection, where the numerically calculated pres- 
sure curves are presented (Figs. 7 and 8). For now let 
us assume that the third term in Eq. ( p6| ) is negligible as 
compared to the first two. Of these two terms, the first, 
P,„, is linear in the density whereas the second term, 
^'hyd, is quadratic. As a result, the relative importance of 
Pm and -Phyd depends on the plate separation d. At large 
d the density in the mid-plane region is small, so that 
Phvd ^ Pm- The main correction to to the PB pressure 
( P7| ) then comes from the change of the mid-plane den- 
sity, c{d/2) — cpB{d/2). Far away from the two plates the 
system behaves as predicted by PB theory with a mod- 
ified, effective surface charge. The mid-plane density is 
depleted relative to PB, since counterions are attracted 
to the vicinity of the charged plates. Hence the pressure 
is smaller than in PB theory. As the plate separation 
decreases and the mid-plane density increases, Phyd can 
become important. 



B. Numerical results 

The general arguments of the previous section can 
be verified by calculating numerically the pressure us- 
ing Eq. ([26[). Fig. 7 shows the pressure as a func- 
tion of the plate separation d for a surface charge 
\a\ = 0.333 C/m^ ~ le/48A^ and bulk ion density 
Cf, = 0.025M (solid line). The pressure is compared with 
PpB (dotted line). The contribution of P,,,,, the first term 
in Eq. (p6|), is also shown (dashed fine). 

The behavior of the pressure at a large range of plate 
separations is shown in Fig. 7a on a semi-logarithmic 
plot. At large d, the pressure is dominated by Pm, as ex- 
pected. It is considerably smaller than the PB pressure, 
due to the reduced effective charge on the plates. At 
lower d the second term in Eq. (p6|), Phyd, becomes dom- 
inant, and the overall interaction is attractive at plate 
separations between 6 and 12 A. Note that the apparent 
sharp decrease in the pressure at a separation of approx- 
imately I3A is artificial, and results from the divergence 
of the logarithmic scale as the pressure approaches zero. 
Fig. 7b shows the same pressure using a linear scale, in 
the region in which it becomes negative (attractive). The 
net pressure crosses smoothly from positive to negative 
values due to a steady increase in the magnitude of the 
(negative) Phyd- At very short separations Pm dominates 
again, and the pressure coincides with the predictions of 
PB theory. 

Fig. 8 shows the effect of the hydration potential for 
a smaller surface charge, \a\ — 0.119C/m^ ~ le/135A . 
In this case and for all surface charge, \a\ < 0.25 C/m^, 
the pressure is repulsive at all plate separations. The 
correction over the PB result is much smaller than in 
Fig. 7, but still significant. At plate separations of ap- 
proximately 5 to 20A Phyd is the dominant contribution 
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to the deviation from PB, and results in a considerably 
reduced pressure. At larger d the pressure is reduced 
mainly because of the change in the mid-plane density. 



C. Comparison with AHNC 

Fig. 9 shows a comparison of the pressure obtaiwjd in 
our model (a) and in the AHNC approximationED (b), 

at surface charge \a\ = 0.333 C/m^ - le/48A^ The 
same short-range hydration potential is used in the two 
calculations. The main plots show the pressure using a 
logarithmic scale. The insets show the pressure on a lin- 
ear scale in the region where it becomes attractive. The 
full pressure (solid line) is compared in Fig. 9a with the 
PB pressure (dotted line). In Fig 9b the AHNC pressure 
(solid line) is compared with the pressure obtained using 
an electrostatic and hard core interaction only (dotted 
line). Since the AHNC approximation accounts for ion- 
ion correlations, there are differences between the pres- 
sure curves in our model as compared to the AHNC 
approximation. However a comparison of Figs. 9a and 
9b shows that very similar qualitative and even semi- 
quantitative effects of the hydration interaction are found 
in the two calculations. 

A comparison for smaller |cr| = 0.119C/m^ ~ 

le/135A is shown in Fig. 10. The solid line is the 
pressure in our model and the dashed line is the AHNC 
pressure. The dotted line shows the PB pressure. As in 
Fig. 9, the qualitative effect is similar in the two calcula- 
tions. 

Since the AHNC approximation takes into account ion- 
ion correlations, the comparison allows us to assess the 
relative importance of correlations and discrete solvent 
effects. The results shown in Figs. 7 and 9 indicate that 
discrete solvent effects can be much larger than correla- 
tion effects induced by the electrostatic interaction. For 
smaller surface charge, as in Figs. 8 and 10, these effects 
are of similar order of magnitude. In the AHNC approx- 
imation the pressure includes an electrostatic term due 
to correlations between ions in the two halves of the sys- 
tem, in addition to the hydration and mid-plane density 
contributions. In Fig. 10 this term is of similar order of 
magnitude as Phyd, and is the main source for the differ- 
ence between the solid line (our model) and dashed line 
(AHNC). For larger surface charge, as in Fig. 9, Phyd 
becomes much larger than the electrostatic contribution. 

When divalent ions are present in the solution, corre- 
laticp. effects become much larger than in the monovalent 
caseEfrE^. Discrete solvent effects are also modified, since 
the effective short-range interaction between two diva- 
lent ions is different in the two cases. When the charge 
on the ions is doubled, the electrostatic interaction be- 
tween two ions increases by a factor of 4. The ion-ion 
separation where the electrostatic interaction is equal to 
ksT increases from Ib — 7 A to almost 30 A. On this 
electrostatic scale, the water molecular size 3 A) is 



much smaller than in the monovalent case. Hence we 
can expect the solvent to be more similar to a continu- 
ous dielectric medium. Indeed, the correction to the X/fr 
potential between two (artificial) Na^+ ions in waterEa is 
found to be purely repulsive, and is significant only at 
separations below ~ 10 A, where the electrostatic inter- 
action is considerably larger than fc^T. Thus we expect 
discrete solvent effects to be less important than correla- 
tions in the divalent case. 



D. Further analysis 

1. Large plate separations 

As discussed above, the hydration term becomes small 
at large d, compared to the change in the mid-plane den- 
sity. In order to study the contribution of the mid-plane 
density to the pressure, let us assume that the plate sep- 
aration d is much larger than all other length scales in 
the system: &, dhydjAu. The two plates are then de- 
coupled and the mid-plane potential can be written as 
^'(d/2) ~ 2*i(d/2) where 5'i(d/2) is the electrostatic 
potential at a distance d/2 from a single plate. We as- 
sume also that ^ h, which is usually the case when 
the surface charge density is large. At a large distance 
from the plate the single plate profile is a PB profile, cor- 



responding to a renormalized surface charge a, 
contrihiition P„ 
foUowsQ: 



offL 



The 

to the pressure can then be written as 



ttIbXJj 



1 - 



26cff 
A 



D 



'd/\D 



(30) 



where bcs — l/27rZs|crcff | is the effective Gouy-Chapman 
length. A similar expression holds for the PB pressure, 
with the nominal Gouy-Chapman length b used instead 
of focff- We thus find that: 



1 - 2bas/\D ^ ^_^bcs-b 



PPB 1 - 2b/ \ 



(31) 



In Ref. |3^ an analytical expression for 6cff — 6 is found. 
Its general behavior is: 



Bt_ 



(32) 



with a numerical prefactor of l/127r in the limit b <C d^yd 
and a numerical prefactor l/47r in the limit b ^ dhyd- 
The parameters of the hydration interaction dhyd — 7 A 

Q 3 c — 

and Bt ~ —500 A are as defined in Sec. [I C . 

A careful treatment of Eq. (p^) shows that the second 



and third terms also add a contribution to the pressure 
that should be regarded as linear in the density, although 
this contribution is small. For large enough d the integra- 
tion range in the second term of Eq. (|2^) can be extended 
to be between —oo and +oo because dB/dz has a finite 
range. In addition all quantities can be replaced by their 
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mid-plane values. We then find that the second and third 
terms of Eq. ( ]26| ) give: 

- fcisT ^ dz' c, (z)c, (z') ^ (z' - z) 

- 2kBTBtcl 

~ ^kBTBtJ2h{d/2) - Cb][cj{d/2) - Cb] 

+2kBTBtCbY,hid/2) - Cb] (33) 

i 

The first term is quadratic in [ci{d/2) — Cb] and can be 
neglected relative to Pm at large d. The second term is 
linear, although small because BjCb ^ 1. It accounts for 
the small difference between the dashed and solid lines 
at large d in Fig. 7a. 

2. Hydration pressure 

The behavior of Phyd, the hydration pressure term, can 
be understood as follows. As a zero-th order approxima- 
tion, the ion density is dominated by electrostatics and 
can be replaced in Eq. ( |29| ) by its PB value. Fig. 11 shows 
that this gives a very good approximation. Hence we can 
write: 

pd/2 pd 

-Phyd ^ / dz / Az' c-p^^i{z)F{z' - z)c-p^,j{z') 

Jo Jd/2 

(34) 

where 

F(z) . -fc.T^ (35) 
dz 

represents the force between two planar ion layers sep- 
arated by a distance z. The following behavior of F{z) 
can be inferred from Fig. 3. At inter-layer separations 
z < dhc = 2.9 A F{z) is positive (repulsive). At larger z 
the value of B{z) increases from its large negative value 
at z = dhc to zero over a few Angstroms, leading to a 
strongly attractive (negative) F{z). A closer inspection 
of Fig. 3 shows that F{z) is oscillatory, due to the local 
maxima and minima of B{z). As we shall see below these 
fine details are smoothed away when two diffusive layers 
of finite thickness interact. 

The behavior of Phyd in Fig. 11 can now be understood 
as follows. Most of the counterions are concentrated 
near the two plates, in layers whose thickness is of order 
b = 1.06 A. Note that b is small compared to dhyd — 7 A. 
When d > dhyd these two layers do not interact directly 
with each other through the short-range interaction. Ions 
in the two sides of the mid-plane interact with each other, 
leading to a negative (attractive) Phyd • As c? is decreased 
towards dhyd, larger and larger ion densities come into 



contact through F(z) and the magnitude of the nega- 
tive -Phyd increases accordingly. The gradual increase in 
the magnitude of Phyd reflects the algebraic decay of the 
density profile near each layer. When d decreases below 
~ 2dhyd — 14 A, the magnitude of Phyd increases more 
rapidly, as the ions in the two layers interact with ions 
in the mid-plane region. 

The behavior of Phyd changes when d decreases be- 
low dhyd- Most of the contribution to Phyd now comes 
from the interaction between the dense counterion lay- 
ers near the two plates. As d decreases these layers are 
separated by correspondingly decreasing distances. The 
hydration pressure follows roughly the structure of F{z). 
It is strongly attractive for d ^ dhc and repulsive for 
d < dhc. The fine details of F{z) are smoothed due to 
the thickness of the diffusive ion layers. 

As the plate separation decreases below dhc towards 
contact Phyd tends to zero, as it should since P(0) — 
—ksT ^|^_Q = 0. One implication of this result is that 
Pm returns to be the dominant contribution to the pres- 
sure, even for high surface charges. Another implication 
is that the short-range interaction becomes unimportant. 
As in PB theory, the ions in the region between the two 
plates become essentially a confined ideal gas, and their 
total number is determined by charge neutrality. Thus 
Pm coincides with the PB pressure matching the nominal 
surface charge density a. This is seen clearly in Fig. 7. 



3. Small plate separations 

In experiments the actual surface charge is usually not 
exactly known, because the number of ions dissociating 
from the surface is uncontrolled. The PB charge is then 
fitted to the large separation behavior. This charge can 
be significantly smaller than the actual surface charge, as 
discussed above. The interpretation of our results is then 
as follows. At plate separations below approximately 2 
nm, an attractive force appears, due to Phyd. This force 
can reduce the net repulsion, or even induce a net attrac- 
tion, depending on the surface charge on the plates. As 
the plate separation decreases below the range of the hy- 
dration interaction dhyd — 7 A, Phyd decreases and even- 
tually tends to zero. The pressure then coincides with 
the PB pressure matching the nominal surface charge. 
As was pointed out in Ref. ^ this leads to an apparent 
strong repulsive force when compared with the PB curve 
fitted to the large separation behavior. As an example, 
the pressure corresponding to ct = 0.25 C/m^ ~ 1 e/64 A 
is shown in Fig. 12 as a function of d (solid line) using 
a linear scale. The dashed line shows the PB pressure 
curve using an effective surface charge chosen to match 
the large d behavior of the solid line. When the two lines 
are compared a strong (apparent) repulsive contribution 
is seen in the solid line below d ~ 5 A, and an attractive 
contribution is seen for 5A< d < 15 A. 
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VI. CONCLUDING REMARKS 

Summarizing our results on the pressure, we find that 
hydration effects can be understood as arising from two 
contributions. The first contribution is the change in the 
mid-plane ion density. This contribution dominates at 
large plate separations and can be understood in terms 
of an effective PB surface charge in our model. The ef- 
fective PB charge is smaller than the nominal charge due 
to the accumulation of counterions in the vicinity of the 
charged plates. Thus the pressure is reduced relative to 
PB theory, using the same surface charge. 

As an alternative viewpoint, the PB surface charge can 
be chosen to match the large plate separation of the pres- 
sure in our model. When this is done, an apparent repul- 
sive force appears in our model at very small plate sepa- 
rations 5 A), as compared with the fitted PB pressure. 

The second contribution to the pressure is the direct 
solvent mediated attraction between ion pairs in the two 
halves of the system. This latter term can become dom- 
inant at plate separations between ~0.5nm and ~2nm. 
It can induce a net attractive interaction between the two 
plates when the surface charge is high. 

Attraction betweeriJ.i|ke-charged surfaces is never pre- 
dicted by PB theoryHa. On the other hand, mecha- 
nisms involving correlations can lead to attraction. Sev- 
eral approaches have shown that ion-ion correlations 
can haiie_t|his effect, in the framework of the primitive 
modelEOtfl. In practice, this attraction can be strong 
enough to overcome the mean field repulsion when diva- 
lent ions are present in the solution. When there are only 
monovalent ions in the solution, ion-ion correlations have 
a much smaller effect. Another mechanism that can lead 
to attraction is the van der Waals force, arising from cor- 
relations between the polarizations on the two surfaces. 
As we find in this work, solvent mediated forces, related 
to ion-solvent correlations, are another mechanism that 
can induce inter-surface attraction. In some cases (mono- 
valent ions, small separation, large surface charge) they 
are the leading mechanism for attraction. 



A stro ps deviation from PB predictions is indeed 
measuredrnrn between charged surfaces in aqueous so- 
lution at separation below ^ 2 nm. The force includes an 
oscillatory contribution, with a period corresponding to 
the water molecular size. This force is due to the struc- 
turing of water in layers between the surfaces. In addi- 
tion to this oscillatory contribution, an additional strong 
contribution isi-seen, which is often referred to as the hy- 
dration forceacll. The aqueous pair potential model of 
Rcf. was a first step towards the understanding of this 
force. A more realistic picture will probably emerge if a 
proper effective ion-surface interaction will be included, 
in addition to the effective ion-ion interaction. In addi- 
tion, the modification to the ion-ion effective potential in 
a confined geometry may also be important. In order to 
assess the importance of these effects, further simulation 
results are needed as an input to the model. 

The aqueous pair potential and the free energy (|^) 
involve various approximations, which are discussed ex- 
tensively in Refs. |3^,^. Nevertheless, the large modi- 
fication to the PB prpsswre^ as obtained also using the 
AHNC approximatioiEilc3't3 , indicates that the solvent 
effects on the ion distributiop-ate a crucial ingredient in 
the origin of hydration forcefS'cJ. The semi-quantitative 
agreement of our results with the AHNC approximation 
indicates that our formalism captures the important ef- 
fects and suggests its further application in non-planar 
geometries, where the AHNC approximation is not ap- 
plicable. 
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APPENDIX A: DERIVATION OF THE PRESSURE 

The free energy of the system is given by the sum: = flpB + Ail with flpB and Afi defined as follows: 

Al] = ifcsTV / dz / dz' c^{z)cj{z')B^j{z' ~ z) (Al) 
^ ,j Jo Jo 

We now imagine that the separation between the two plates is increased from d to d + Sz by adding a 'slice' of width 
Sz between the planes zq and zq + 6z. We map the regions < z < zq and zq < z < d in the original system to the 
regions < z < zq and zq + Sz < z < d + 5z in the modified system, respectively. We then have: 
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C 



(A2) 



The first term can be integrated by parts. With the conventions described above, the boundary terms can be 
written as foUows: 



^6 



d* 



, , d* 





dz 



—5 [zq) - -^-(^o) -. (20 + oz) + --7-(za) 

dz dz dz dz 

d^* 

= -'^{zo)^izo)6z 



(A3) 



where use of the boundary conditions at z = and z — d has been made. Using this relation and the Poisson equation 
we obtain: 



i 

■Jo 



e fd-^y 



TT V dz / 



ksT c, I In - 1 



(A4) 



To compute (5Ai7, Ai7 can be separated to the following three terms 



1 rzo rzo 

Arj = -fci3rV/ dz / dz'c^(z)cj(z')By (z' - z) 

2 . ^. Jo Jo 

+ JfcBr^ [ dz [ dz' c,{z)cj{z')B,Az' - z) 

^ , , Jzo J Zo 

d 



rzo ra 

+kBTY^ dz dz' c^{z)cJ{z')BiJ{z' - z) 

, , Jo Jzn 



The variation of Ci in these three terms gives: 



pd pd 

SAni = kBT'Y dz dz' c,(z)Scj{z')B,j{z' 



z) 



(A5) 



(A6) 



The variation of the third term in Eq. ( |A5| ) gives two additional contributions, one from the variation of B{z' — z) 
under the insertion of the 'slice' at zq: 



SAn2 = Sz-kBTy2 dz / dz'c.(z)cj(z')-r^(z'-z) 
„ "'0 Jzo dz 



(A7) 



and the other from the integration over the 'slice' itself: 



Pd 

^ Sz ■ ksTy^ / dzCi(zo)cj(z)i3y (z - Zo) 
. Jq 



(A8) 



Summing up all the contributions to SCI we have: 
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sn 



6Ani + SAa 




2 

) (zo) + e,*(zo) + fcsT 



dz 



-bz kgT 



' Sciiz) X { ej*(z) + ksTln 



Ciiz) 



dz 



dz 



Using the equilibrium equation this reduces to: 

]ci{zo) - 



sn , 

— = kBT^ 

oz 



e 

8^ 



dz 



In 



ksTY^f 



dzcj{z)Bij{z - zo) 



k„T 



j 



dz' Cj{z')B,j{z' - z) 



dz' c,iz)cj{z')^iz' ~ z) 



dz' Ci{z)cj{z' 



^-r^{z -z) 
dz 



This result can be readily generalized to the case of several ion species, as in Eq. 



(A9) 
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Fig. 1: Short-range effective potential between Na ion pairs, adaated from Ref. |27| using simulations in a bulk 
NaCl aqueous solution of concentration 0.55M, at room temperatureES. The potential is shown in units of ksT, as a 
function of the distance between the ion centers. The Coulomb interaction is subtracted to show only the short-range 
hydration effect due to the water molecules. For ion-ion separations below 2.9A a hard core interaction is taken. 

Fig. 2: Schematic description of the pair potential model. An aqueous ionic solution confined between two charged 
plates in (a) is replaced by ions in a continuum dielectric medium with electrostatic and short-range interactions 
Uij{r) = Uij{\r\) in (b). The coordinates z = and z = d designate the contact positions of the ions with the plates. 
The distance of closest approach is equal to dhc/2, where dhc is the hard-core diameter of the ions. 

Fig. 3: The effective layer-layer interaction B{z) in a planar geometry, as obtained from the potential of Fig. 1 using 
Eq. (H) (solid line). The oscillating structure of the radial potential shown in Fig. 1 is apparent in the secondary 
minima of B{z). 

Fig. 4: Counterion density profile (solid line) obtained from numerical solution of Eq. ^ with the hydration interaction 
as of Fig. 3, plotted on a semi-log plot. The bulk ion density is Cb = 0.025M and the surface charge is |cr| = 0.333C/m^. 
The dielectric constant is £ = 78 and the temperature is 298K. The distance between the plates is d = 50A. The 
density profile is symmetric about the mid-plane at z = 25 A. The dotted line shows the corresponding density profile 
obtained from the PB equation. The symbols ('x') show the density profile obtained in the AHNC approximation, 
using the same parameters. 

Fig. 5: The ratio of the positive ion density obtained from Eq. (|^) and the value obtained from PB theory, for surface 
charges \cr\ — 0.333C/m^ (dashed line), O.lC/m^ (solid line) and 0.0333C/m^ (dotted line). All other parameters are 
as in Fig. 4. 

Fig. 6: The ratio between the positive ion density obtained from Eq. and its PB value, for plate separations d 
equal to (a) 50 A (solid line), 35 A (dashed line), 20 A (solid line), (b) 10 A and (c) 5 A. All other parameters are as 
in Fig. 4. Each curve is shown between the plate at z = and the mid-plane z = d/2. 

Fig. 7: (a) Pressure between two plates with surface charge \a\ = 0.333C/m^, as a function of the plate separation 
d, on a semi-log plot. All the parameters are as in Fig. 4. The solid line shows the overall pressure P obtained from 
Eq. (^^ . The dashed line shows the contribution Pm resulting from the mid-plane density and the dotted line shows 
the PB pressure, (b) The same curves on a linear scale, in the region where the overall pressure becomes negative, 
i.e., attractive. 

Fig. 8: The repulsive pressure between two plates with surface charge \a\ = 0.119C/m^, as a function of the plate 
separation d. All other parameters are as in Fig. 4. The solid line shows the overall pressure P, the dashed line shows 
the contribution Pm of the mid-plane density, and the dotted line shows the results of PB theory. 

Fig. 9: Comparison between the pressure obtained (a) in our model and (b) in the AHNC approximation, using 
the same short-range hydration potential (solid lines). All the parameters are as in Fig. 7 (|ct| = 0.333C/m^). The 
pressure is shown as a function of the plate separation d. A semi-logarithmic scale is used in the main plots and a 
linear scale is used in the insets. In (a) the dotted line shows the PB pressure. In (b) the dotted line shows the pressure 
obtained in the AHNC approximation when the ion-ion interaction includes only the hard core and the electrostatic 
interactions. 

Fig. 10: Comparison between the pressure obtained in our model (solid line) and the AHNC approximation (dashed 
line), for a surface charge jfjl = 0.119C/m^. All the parameters are as in Fig. 8. The pressure is shown as a function 
of the inter-plate separation d using a semi-logarithmic plot. The dotted line shows the PB pressure. 

Fig. 11: The hydration pressure Phyd as a function of the plate separation d (solid line). All the parameters are as 
in Fig. 4. The dashed line shows the approximation to Phyd obtained by replacing the ion density in the integral of 
Eq. ( p9| ) by the PB ion density. 

Fig. 12: Repulsive pressure between two plates with surface charge \a\ = 0.25 C/m^, as a function of the plate 
separation d, using a linear plot (solid line). All other parameters are as in Fig. 4. The pressure is compared with the 
PB pressure curve fitted to the large separation behavior, with \(Tcs\ — 0.09 C/m^ (dashed line). 
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